GLOBAL ASYMPTOTIC SOLUTIONS FOR 
NON-RELATIVISTIC MHD JETS AND WINDS 



Jean Heyvaerts 
Universite Louis Pasteur, Observatoire de Strasbourg 1,4 
^ . heyvaert@astro.u-strasbg.fr 

o 

C\| 1 and 

CD ' Colin Norman 

C/3 ■ 



> 
00 
(N 



Department of Physics and Astronomy, Johns Hopkins University 
and Space Telescope Science Institute 2,2, 



normanOstsci . edu 



On 

O ! ABSTRACT 

m 
o 

^ ' We present general and global analytical solutions, valid from pole to equator, 

for the asymptotic structure of non-relativistic, rotating, stationary, axisymmet- 
ric, polytropic, unconfined, perfect MHD winds. The standard five lagrangian 
first integrals along field lines are assumed known. 

The asymptotic structure of such winds consists of field-regions virtually de- 
void of poloidal current. We show that an Hamilton- Jacobi equation, or equiv- 
alently a Grad-Shafranov equation, gives the asymptotic structure in the field 
regions. These field regions are bordered by current-carrying boundary layers 
around the polar axis and near null magnetic surfaces. Current closure is achieved 
in a number of separate cells bordered by null surfaces. The solution is given in 
the form of matched asymptotics separately valid outside and inside these bound- 
ary layers. The polar boundary layer is pressure-supported against the pinching 
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force exerted by the axial poloidal current and has the structure of a current 
pinch, while the null-surface boundary layers have the structure of current sheet 
pinches. We establish a consistency relation between the residual poloidal cur- 
rent at large distances and the axial pressure. We find a similar relation for the 
current sheets at null surfaces. We further consider the case where the polar 
boundary layer is force free. The geometry of magnetic surfaces in all parts of 
the asymptotic domain is explicitly deduced in terms of the first-integrals. 
The solutions have the following general properties: 

I. For winds which are kinetic-energy dominated at infinity we derive WKBJ 
analytic solutions whose magnetic surfaces focus into paraboloids. The current 
slowly weakens as the inverse of the logarithm of the distance to the wind source 
while the axial plasma density falls-off as a negative power of this logarithm. 

II. For winds carrying Poynting flux at large distances the solutions asymp- 
totically approach to nested cylindrical and conical magnetic surfaces. 

Subject headings: jets,MHD 

1. Introduction 

We have previously established (Heyvaerts & Norman 1989) that any stationary ax- 
isymmetric magnetized wind will collimate at large distances from the source, under perfect 
MHD conditions and polytropic thermodynamics, to paraboloids or cylinders along the sym- 
metry axis according to whether the electric current (or Poynting flux) brought by the wind 
to infinity asymptotically vanishes or is finite. Our result however did not discriminate be- 
tween these two possibilities, nor described any global asymptotic solution. The aim of this 
paper is to partially fill this gap. We consider here non-relativistic winds and concentrate on 
describing analytically their structure in the aymptotic region, for a supposedly given set of 
first integrals and for different, a priori possible, values of the circumpolar current at large 
distances. 

We show here that the asymptotic structure of rotating MHD winds consists of vast re- 
gions where the poloidal current density is negligibly small, bounded by thin regions where 
residual asymptotic poloidal current flows. These regions have at large distances the char- 
acter of boundary layers. They are located in the vicinity of the polar axis and of the null 
magnetic surfaces. We obtain solutions valid in all these regions separately and produce a 
global solution by asymptotic matching. 

The specific topic which we are pursuing here is the construction of a general, non-self- 
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similar, asymptotic solution globally valid from the polar axis to the equator, for a given set 
of supposedly known first integrals. 

There has been in the past few years considerable progress in this field, both in the 
derivation of special exact solutions to the wind equation and in numerical solutions to 
them, both time dependent and stationary. Focusing of magnetized winds appears to be 
a robust property of rotating MHD winds (Blandford & Payne 1982; Heyvaerts & Norman 
1989). Most analytical solutions involve some sort of self-similarity (Lovelace et al. 1991; 
Contopoulos & Lovelace 1994; Ostriker 1997; Trussoni et al. 1997; Vlahakis & Tsinganos 
1998, 1999). Lynden-Bell (1996) has constructed quasi-static, force-free collimated structures 
that arise naturally from a wound up magnetic field pushing out from a disk. The dynamics 
of such structures has been studied by Kudoh et al. (2002). 

Sauty et al. (1992a,b), Sauty & Tsinganos (1994) and Sauty et al. (1999) have analyt- 
ically considered particular models of non-polytropic winds and found that non-cylindrical 
asymptotics can be achieved only when magnetic pinching is negligible and there is over- 
pressure in the vicinity of the axis. 

Shu and collaborators have extensively developed a X-wind model for outflows and in 
Shu et al. (1995) (plus references therein) they actually solved a Grad-Shafranov equation for 
the shape of the magnetic surfaces (Shafranov 1996). They have an inner boundary condition 
of finite pressure. The launching region for T-Tauri winds has been examined recently by 
Anderson et al. (2003). Interesting disk instability can be initiated by wind driven angular 
momentum loss Cao & Spruit (2002). Tomisaka (2002) has found that the collapse of rotating 
magnetized molecular clouds and the resulting bipolar outflows are inextricably related. 

Transfield force balance has also been studied by numerical simulations. These simu- 
lations studied hoop stresses collimation in different parameter regimes. Ustyugova et al. 
(1995) and Romanova et al. (1997) obtained two-dimensional solutions for jets emitted by a 
source with Keplerian rotation and confirmed the focusing effect of the hoop stress, as did 
Ouyed and Pudritz (1997a); Ouyed & Pudritz (1997b, 1999) who also found time-dependent 
behaviour. Lery & Frank (2000) found by two-dimensional simulation that winds orginating 
from a disk with a Keplerian rotation profile have a dense, current-carrying, central core 
surrounded by an almost current-free region. Ustyugova et al. (1999) numerically calcu- 
lated this same problem using a model with a hot wind source region in the vicinity of the 
polar axis and found very little collimation within the simulation region. Ustyugova et al. 
(2000) studied formation of collimated Poynting jets associated with an uncollimated hy- 
dromagnetic outflow. Bogovalov & Tsinganos (1999) numerically found collimation to be 
most effective for a particular class of objects which they describe as "efficient rotators". In 
a subsequent paper Tsinganos & Bogovalov (2000) discussed the efficiency of collimation in 
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the case of the solar wind. Boundary conditions at the base of the flow were found to be 
important. Bogovalov (1996, 2000) studied cold non- stationary flow from a monopole and 
found significant Poynting flux conversion to kinetic energy. 

Since Heyvaerts & Norman (1989), papers similar in spirit to our own research have 
appeared. Li (1993) and Begelman & Li (1994) studied kinetic energy dominated winds and 
obtained asymptotic solutions that agreed with the paraboloidal solutions we had previously 
obtained. Tomimatsu (1995) constructed solutions in different regions including the pole and 
the equatorial null surface quite similar in spirit to the analysis presented here. Okamoto 
and collaborators (Okamoto 1999, 2000; Beskin & Okamoto 2000) emphasized the issue of 
current closure and its effect on the geometry of the solution. These issues which are, in 
fact, important only in boundary layers such as the equatorial current sheet, are now fully 
treated in this paper. 

Our presentation is structured as follows. Section (2) reviews the basics of rotating 
MHD winds. Section (3) deals with the field-regions where almost no poloidal current flows. 
For Poynting jets the transfield equation becomes an Hamilton-Jacobi equation (Eq.(26)), 
or equivalently a Grad-Shafranov equation (Eq.(18)), which we solve in section (4) (Eq.(30)). 
Kinetic winds are solved using a WKBJ approximation. In section (5) we obtain the solution 
for the polar boundary layer (Eqs.(43)-(44)). This solution, which is similar to a Bennet 
pinch, is then matched to the field-region solution and a relation between the asymptotic 
current and the axial pressure is derived (Eq.(50)). The case of the force-free polar boundary 
layer is also discussed, giving a relation similar to the standard Bennet pinch (Eq.(70)). The 
matching procedure specifies the density along the polar axis (Eqs.(54)-(50)), the asymptotic 
circumpolar current and the radius of the current-carrying region (Eq.(57)). A slow loga- 
rithmic decline of the axial density and current is found in case of kinetic-energy-dominated 
winds, justifying the WKBJ approach. In section (6) we similarly obtain a solution valid 
in the vicinity of a null magnetic surface (Eqs.(76)-(79)) and match to the field-region so- 
lution. The gas pressure at the null surface balances the toroidal magnetic pressure outside 
the sheet (Eq.(102)). In section (7) the shape of the magnetic surfaces has been calculated 
in all regions and in both regimes. In the paraxial field-region of kinetic winds it is given 
by Eqs.(95)-(94). In the polar boundary layer it is given by Eq.(100). Inside an equatorial 
null-surface boundary layer of a kinetic wind it is given by Eq.(107), while inside the equa- 
torial null-surface boundary layer of a Poynting jet it is given by Eq.(106). Our conclusions 
regarding the general properties of non-relativistic rotating MHD winds are presented in 
section (8). 
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2. Axisymmetric Stationary MHD Flows 

2.1. Notation and Definitions. 

The formulation of stationary axisymmetric rotating MHD winds have been presented 
in a number of papers (Weber & Davis 1967; Okamoto 1975; Mestel 1999; Sakurai 1985; 
Heyvaerts & Norman 1989; Heyvaerts 1996). Cylindrical coordinates r, 9, z are used with 
unit vectors e r , e$, e z . The spherical distance of a point to the origin is denoted by R. Vectors 
are split into toroidal and poloidal parts, indicated by subscripts 9 and P respectively. The 
MKSA system of units is used with p being the magnetic permeability of free space. The 
poloidal magnetic field can be expressed in terms of a magnetic flux function a(r,z), such 
that 

-> 1 da _ 1 da _ 

B P = — Tre r + -—e z . (1) 
r oz r or 

The magnetic flux $ m through a circle of radius r centered on the axis at altitude z is 
$ m = 2na(r,z). The magnetic surfaces, generated by the rotation of field lines about the 
axis, are surfaces of constant value of a(r, z). The value of a labels magnetic surfaces or field 
lines. 



2.2. First Integrals. 

Mass conservation implies 

pv P = a(a)B P , (2) 
where a (a) is a first integral. Flux freezing implies that 

pv e = a(a)B e + prVL(a), (3) 

where Q(a) is a second first integral. Conservation of specific angular momentum implies, 

Conservation of the specific energy E(a) implies: 

rn/\ 1 2 1 2 IP*. rB e Q(a) 

E(a) = -v 2 P + -v 2 e + -!— y - + $ G tAJ. 5 

2 2 7 — 1 p p a{a) 

E(a) is the total specific energy carried by the wind on magnetic surface a, in kinetic and 
Poynting flux form. $ G is the gravitational potential. Equation (5) is the Bernoulli equation. 
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The thermodynamics is described by a lagrangian polytropic law, relating pressure and 
density, of the form 



This defines a fifth first integral, Q(a), the polytropic entropy. Equation (6) may represent 
adiabatic or more complex thermodynamics. 



We define terms that will be used in this paper. The asymptotic domain consist of all 
points which are located on their own magnetic surface far away from the Alfven point. A 
magnetic surface is said to be asymptotically parabolic if r and (z/r) both become infinite. 
This does not imply that z should vary as a power law of r at fixed a. A magnetic surface 
is said to be asymptotically conical if r becomes infinite at large distances while (z/r) ap- 
proaches a finite limit tan^^a). This does not imply that (z — rtanip^) should approach 
a finite limit. Conical magnetic surfaces may have parabolic branches. An asymptotically 
cylindrical magnetic surface is one on which the axial distance r approaches a finite limit 
roo(a). 

When the asymptotic magnetic structure consists of cylindrical surfaces nested inside 
conical ones, the cylindrically focused one are refered to as forming the jet, and the other 
surfaces as forming the conical wind. The jet edge is defined as being the magnetic surface 
which separates asymptotically cylindrical from asymptotically conical surfaces. This termi- 
nology does not imply that only cylindrical asymptotics should give rise to structures that 
would appear to the observer as an astrophysical jet. 

A neutral or null magnetic surface is one along which the poloidal field vanishes. The 
toroidal field then also vanishes on it. The immediate vicinity of neutral magnetic surfaces is 
of special interest, being regions where electric currents flow in the asymptotic domain. We 
refer to these regions as neutral surface boundary layers. Similarly, the vicinity of the polar 
axis is a region of electric current flow. We refer to it as the polar boundary layer. Outside 
the boundary layers is the field-region. 

The space between two neighbouring neutral magnetic surfaces is refered to as a cell. 
The total electric current through a circle of axis z passing through the point (r, z) is: 



v = Q(°)p 



(6) 



2.3. 



Semantics. 




(7) 



J vanishes at neutral surfaces and is negative for positive Q and a. Therefore we use the 
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quantity 



which is positive for positive f2 and a and proportional to the current J. We refer to it as 
the current, noting that the true physical current is J = —2ttI. 

Strucures which convey a finite circumpolar electric current to infinity are refered to as 
Poynting jets. Kinetic winds carry no Poynting flux, or current, at infinity. 



2.4. Bernoulli Equation 

The toroidal components of the velocity and of the magnetic field can be obtained from 
the angular momentum conservation and isorotation law in terms of p and r as 

L p L-r 2 Q 
v 9 = - + - 5 , 9 

r r p ct z — p 

B e = ^ — . (10) 

r p a 2 - p 

For well-behaved solutions the numerators must vanish when the denominators vanish. The 
cylindrical radius r equals the Alfven radius, r^a), defined by Q(a) r 2 A (a) = L(a), when the 
density p equals the Alfven density, Pa(cl) = Poa 2 {a). The Alfven radius defines the position 
of the Alfven point on field line a. Asymptotically, for r ^> r^, the azimuthal velocity 
vanishes while the variable / = —rBe/po approaches the value 

i (ii) 

p a(a) 



2.5. Transfield Equation 

The projection of the equation of motion on the normal to the magnetic surfaces gives 
the transfield equation, or generalized Grad-Shafranov equation (eq. (A-l) of appendix A). 
It determines the shape of magnetic surfaces in terms of the first integrals and of the density. 
The set of coupled Bernoulli and transfield equations is to be solved. As shown in appendix 
A, the transfield equation can be reexpressed such as to give the curvature (dip/ds) of poloidal 
field lines in terms of various surface functions and of the current variable /. The resulting 
form of the transfield equation is: 

( , P \ 2 dip 1 Va - / |Va| 2 ^ \ Va - 
\ pa) ds p|Va| \2/i r^ J \\/a\ 
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- (L - I /a? 1^ + J^- • V («£) (12) 
v 7 ; r 3 |Va| pr 2 |Va| \ 2 / 

The terms on the right are the components of the poloidal magnetic pressure gradient, gas 
pressure gradient, gravity, centrifugal force and hoop stress respectively perpendicular to 
magnetic surfaces. The hoop stress is proportional to jp x Bg and therefore vanishes with /. 



3. Field Regions 

3.1. Asymptotically Dominant Forces. 

In the large- ,2 limit, Eq.(12) simplifies. The poloidal velocity \vp\ reaches a terminal 
value and the curvature term on the left of equation (12) approaches zero. A definite ordering 
between the terms on the right of equation (12) occurs on magnetic surfaces on which both 
z and r tend to infinity. On asymptotically cylindrical magnetic surfaces, r approaches a 
finite limit r 00(a). However, when r 00(a) ^> r^a) it still remains possible to use the large-r 
approximation. In the large r, large z and small (p/ Pa) limit, noting that pr 2 is bounded 
from above (Heyvaerts & Norman 1989) and that the gravity term declines as 1/r 2 while the 
centrifugal force term declines as 1/r 3 , Eq.(12) reduces to: 



ds p al |Va| \2p r 2 J a |Va 

— * 

When I approaches a finite value, the hoop stress term, proportional to ft ■ VI decreases 
on conical magnetic surfaces as (1/r) since n ■ W « I jr. The hoop stress force then 
dominates all the other terms on the right of (13). The poloidal field Bp decreases in this 
case as (1/r 2 ). The gradient of the poloidal magnetic pressure then drops with axial distance 
as (1/r 5 ), so that the corresponding term in Eq.(13) declines as (1/r 3 ). According to Eq.(2), 
the density decreases as Bp when the velocity has reached its terminal value so that the gas 
pressure gradient scales as (l/r 27+1 ) and the pressure term in Eq.(13) declines as (1/r 27-1 ). 
This decline is slower than that of the poloidal magnetic pressure if the polytropic exponent 
7 is less than 2, which we assume. The poloidal magnetic pressure is negligible compared to 
the gas pressure except in the ultra-cold limit where Q vanishes. Therefore, we retain the 
gas pressure on the right of Eq.(13). The case of cold winds in which the poloidal magnetic 
pressure is retained is discussed in section (5.5). 

Kinetic winds are characterized by r|Va| and pr 2 approaching zero. The hoop stress 
scales as V/ 2 ~ p 2 r 3 . The gas pressure gradient term scales as p 7 /r, and the poloidal 
magnetic pressure force ~ p 2 jr. Again, gas pressure dominates over poloidal magnetic 
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pressure, except for very cold winds. The hoop stress term dominates over the gas pressure 
force if p 2 r A p 7 . This is so in the kinetic wind solutions derived below, since pr 2 is found 
to decline only as an inverse power of (lnr). 



3.2. Transfield Equation 

Neglecting pressure, the transfield equation (13) then further reduces to: 

2 dip Vl Va -± T A , 
v 2 p~r = +- pz— -V/ 14 

as a |Va| 

Eq.(14) is equivalent to Okamoto's conclusion that j\\B t = pv 2 /R c , R c being the curvature 
radius of poloidal field lines (Okamoto 1999). Note however that the centrifugal force asso- 
ciated with the poloidal field line curvature (on the left of equation (14)) is also negligible. 
Since the poloidal current is bounded from above, the right of equation (14) could a priori 
scale as 1/r, but the curvature (dip/ds) must decline with r faster than 1/r. Otherwise 
the poloidal field lines would not be well-behaved (Heyvaerts & Norman 1989). The left of 
Eq.(14) is therefore of order (v 2 P /R c ) where R c is the curvature radius, much larger than r. 
To lowest order in the small parameter (r/R c ) Eq.(14) therefore reduces to the vanishing of 
its right-hand-side and becomes: 

Va-W = (15) 

This statement has been presented in Heyvaerts & Norman (1989) as a solvability condition 
at infinity. It implies that almost no poloidal electric current density is present in these 
regions. 



3.3. Existence of Current Carrying Boundary Layers and Electric Circuit 

Any physical current system must be closed. Current closure causes I to depend on a 
wherever current is flowing. This implies that some terms of Eq.(13) should balance hoop 
stress where current flows. Since this is not possible wherever V ~ 1/r, this implies that 
regions of current flow must take the form of boundary layers in which gradients are large. 
Note that the hoop stress vanishes where the current vanishes. Boundary layers must then 
occur in the vicinity of the polar axis, where Bg vanishes for finite current density at the 
axis, and near neutral surfaces, where the hoop stress vanishes with B e . Indeed, the toroidal 
field component is generated from the poloidal one by the rotation, but the poloidal field 
vanishes at a neutral surface. In the case of a dipolar symmetry the equatorial plane is a 
neutral surface. 
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We then conclude that Eq.(15) applies in the field-region of the asymptotic domain, ex- 
cept in the boundary layers. Eq.(15) integrates along any trajectory orthogonal to magnetic 
surfaces as long as no null surface (in the vicinity of which equation (15) ceases to be valid) 
is met. The result is: 

I = Ioo(b) (16) 

where b is a label for an orthogonal trajectory. To be specific, we define it as being the value 
of the z-coordinate on this orthogonal trajectory at the polar axis. An extra label, which we 
omit, should precise on which segment of the orthogonal trajectory Eq.(16) applies. Note 
that the magnetic surface label, a, and the orthogonal trajectory label, b, could constitute 
a set of orthogonal coordinates in the poloidal plane that could be substituted to r and z, 
with a playing the role of an angular variable and b the role of a radial distance variable. 

Eq.(15) states that the component of the poloidal current density parallel to Bp vanishes. 
Since 1 00(b) approaches a constant (possibly zero), so does also the other component of jp. 

— ♦ — * — * — * -* 

Thus, the hoop stress j P x B , approaches zero, as does j g x B P , since \B P \ declines rapidly 

— * — * 

with distance. Hence, the Lorentz force j x B approaches zero. This does not imply that the 
field becomes force free, because this decline is due to both current and field approaching 
zero. It is in fact quite interesting to analyze how the different components of the Lorentz 
force individually approach zero at large distance. This requires, however, a more precise 
knowledge of the solution than we have at this point. Therefore this discussion is postponed 
to appendix J. We emphasize that the asymptotic regime discussed in this paper is not 
force-free. 

According to Eq.(16), the enclosed poloidal current is approximately constant through 
any circuit drawn on surfaces orthogonal to magnetic surfaces and running between two 
successive boundary layers. The current enclosed in a circle of increasing radius drawn on 
one such surface would increase from zero on the polar axis to some value at the outer 
edge of a polar boundary layer. The current remains almost constant up to the next null- 
surface, at which it returns to zero (Fig 1) through a boundary layer, vanishing at the null 
surface. Needless to say, Eqs.(15) and (16) are only approximate results. We do not mean 
that the current density is strictly zero in field regions, but only that it is so small that the 
total electric current circulating between two successive boundary layers through a surface 
of constant-6 is much less than the current in the boundary layers themselves. Nevertheless 
a small current in field regions must be present for loo to be a (slowly varying) function of 
b. This picture generalizes for any number of null surfaces and the current system consists 
of more than just two cells (Fig 2), in each of which the electric circuit separately closes. 
1 00(b) approaching zero as b — > 00 means that the poloidal current circuit closes in each 
cell at a finite distance. By contrast, 1 00(b) approaching a finite value, loo, implies that this 
circuit would close at infinity. In reality, the wind has not been blowing for an infinite time 
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and it is externally bordered by a time- dependent region and an expanding shock system. 
This region plays the role of the "region at infinity" in the present stationary model: this is 
where the residual electric currents flow from the circumpolar region to the ret urn- current 
boundary layers shown in Fig. (2). 

Regions about the polar axis and the null surfaces must have the character of thin 
boundary layers since in the asymptotic domain, the gas pressure is a subdominant force 
that can balance Lorentz forces only in a small region about those special places where 
the latter vanishes. These boundary layer sheets must be thin and there must be toroidal 
magnetic pressure equilibrium accross them. The total poloidal current then only changes 
sign at their crossing as Be reverses. 

Our aim now is to solve the wind equations in the asymptotic domain, both in field- 
regions and in current-carrying boundary layers and obtain the resulting shape of the mag- 
netic surfaces. 



3.4. Asymptotic Grad-Shafranov Equation 

From Eq.(14), it is possible to restate Eq.(15) as an equation similar to the familiar 
Grad-Shafranoff equation of magnetohydrostatics. Using the identity (A-2) of appendix A 
and Eqs.(ll), (1) and (2), the equation (14) becomes: 

The above discussion indicates that the curvature radius should be much larger than r, 
so that the left of Eq.(17) is negligible. Expanding the divergence term and using Eqs.(2), 
(1) and (11) to express r|Va|, Eq.(17) is eventually brought to the form of a quasi-linear 
elliptic equation: 

r^=l(^Qf] (18) 




da \ 2tt 2 

The boundary conditions to Eq.(18) are that a = along the polar axis and that a = A on the 
equatorial plane. These are consistent, when / is constant and non-zero, with a depending 
only on the latitude angle ip. In this case Eq.(18) becomes an ordinary differential equation 
for a(ip), which reduces to the form of Eq.(33). Eq.(18) is equivalent to the other forms 
obtained above, in particular it implies Eq.(15). It is also equivalent to the alternative forms 
derived below, in particular Eq.(26). Using Eqs.(ll), (1) and (2), an equivalent form of 
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Eq.(18) is in fact found to be: 

Aa = Va- V (ln(V|Va|)) (19) 
which can be brought by denoting the normal unit vector to poloidal field lines by n = 

— * — * 

— Va/|Va|, to the form 

div = (20) 



4. Solutions in Field Regions 

4.1. An Hamilton- Jacobi Equation 

Noting that the velocity on each surface reaches a terminal velocity foo(o), Eqs. (2), (5) 
and (10) reduce in the asymptotic limit to: 

prVoo(a) = a|Va| (21) 

/ = +^ (22) 
p a 

E = ^k + — (23) 
Eliminating and p between Eqs. (21), (22) and (23), an expression of r|Va| is obtained: 



r |Va| = ^ — (24) 

This relation can be integrated following orthogonal trajectories to magnetic surfaces. The 
curvilinear abcissa along them is denoted a and conventionally increases from pole to equator, 
so that |Va| = da /da. Eq.( 24) becomes 

da Q(a)da . . 

= (25) 



r fi I(a,b)V2^E{ a) — I(a, b)Q(a)/ 



a\a) 



In field-regions, Eq. (25) simplifies, since I (a, b) becomes a function Ioo(b) independent of a. 
It can be restated as: 

lT7 | _ //o/ 00 (6)v / V£(a) - ggjgRTgg (0R . 

r|Va| " n(o) { ' 

If Joo(^) were to approach a non- vanishing value 1^ at large distances from the wind source, 
this equation would explicitly give the flux distribution in space in the asymptotic domain. 
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Particular versions of Eq.(26) for cylindrical or conical magnetic surfaces have in fact been 
obtained and solved in our earlier work (Heyvaerts & Norman 1989). Eq. (26) improves on 
this by not being restricted to either one of these specific geometries. It implies no a priori 
constraints on the structure of the solution and keeps a fully 2-dimensional character. In the 
asymptotic field-region, far from the neutral surfaces and from the polar axis, the transfield 
equation reduces to Eq.(26) with constant Ioo(b), which is of the Hamilton- Jacobi type: 

r\Va\ = /(a), (27) 

By defining S(a) = f^da'/f(a'), it can be converted into 

IVS1 = -, (28) 
r 

— * — * 

This equation can be restated as VS* = —n/r, implying that V x (n/r) = 0. This, with 
Eq.(20), indicates that S should be an harmonic function. The solution for S represented 
by Eq.(33) is harmonic. 

If the poloidal electric current Ioo(b) declines towards zero at large distances, the function 
on the right of Eq.(27) also depends on b. Eq.(26) then does not provide the value of the 
modulus |Va| as a function of a and r (as does Eq.(27)) because it is not known how Ioo(b) 
declines with b. This difficulty can however be circumvented if this variation is so slow that 
it can be treated, as we do below, by a WKBJ type of approach. Note that the case when 
-foo(fr) approaches a non-zero limit, can be dealt with by a WKBJ procedure as well, since 
/oo(6) is not strictly constant, but slowly evolves towards its limit. Important features of 
the asymptotic structure, in particular the difficult question of the connexion between the 
asymptotically cylindrical and the asymptotically conical regions, must be approached by 
considering the non-constancy of Ioo(b) as it approaches its limit. 

Particular solutions to the Hamilton- Jacobi equation (27) for constant Ioo(b) can be 
found by the method of separation of variables. These are described in appendix D, although 
the boundary conditions associated to some of them are different from those we are facing. 
These solutions however exemplify a number of typical structures to which Eq.(27) may 
give rise. General solutions associated with the boundary conditions relevant to the case of 
astrophysical MHD winds can be found by reduction of Eq.(28) to a ray-tracing problem. 

4.2. Ray- Tracing 

Equation (28) is of the form \VS\ = N(f), which is the Eikonal equation for the prop- 
agation of waves in a medium of index of refraction N(r). The wave fronts are represented 
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by surfaces of constant S(r). They may be found by ray-tracing methods, the rays being the 
orthogonal trajectories to the surfaces. So, finding the general solution of Eq.(28) is equiva- 
lent to solving Snell's refraction equation. We show in appendix (B) the equivalence of the 
Hamilton- Jacobi equation with Snell's equation. The polar and the equatorial lines, being 
field lines, must be lines of constant S. The boundary condition to Eq.(28) is that rays start 
perpendicular to the pole and the equator. In the vicinity of a null surface (e.g. the equator), 
the simplified asymptotic transfield equation (15) is invalid. We show in section (7.5) and in 
appendix I that the angle of the magnetic field lines to the equator at the outskirts of this 
boundary layer decreases to zero with increasing b, justifying the view that rays from the 
field- region must nevertheless end perpendicular to it. Similar considerations apply to the 
vicinity of the pole and to other null surfaces: rays must also cross them perpendicularly on 
both sides. 

This optical analogy can be used to find the general solution Eq.(28) by writing the 
appropriate form of Snell's law. The gradient of the refraction index is in this case perpen- 
dicular to the polar axis. Let % be the angle of incidence between this radial direction and 
the tangent to the ray. The ray-tracing equations can be written as: 

Sill % 

dr = —da cos? dz = da sini = k (29) 

r 

where A; is a constant. Unlike the boundary condition at the equator, the condition that i 
approaches zero at the polar axis is not restrictive because N(r), being equal to 1/r, diverges 
there. Equations (29) can be integrated to: 

(z -h) 2 + r 2 = ± (30) 

where h and k are two integration constants. Eq.(30) represents a circle centered at z = h 
on the axis, with a radius D = 1/k. Lines of constant S associated with the solution of 
equation (28) are orthogonal trajectories to a collection of circles centered on the polar axis. 
The boundary condition that the rays connect perpendicularly to the equator implies that 
(h/D) should asymptotically vanish. It is important to note that this is not an exact, but only 
an approximate, asymptotic result. This is consistent with our earlier results (Heyvaerts & 
Norman 1989). We show in appendix C that other conceivable current-enclosing geometries 
are not consistent with a source subtending a finite flux. We have assumed above that a 
unique cell extends from pole to equator. The extension to winds with a larger number of 
null surfaces is described in appendix E. 
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4.3. WKBJ Approximation 

For Ioo(b) approaching a non-zero value as b grows, the field- region of the flow consists 
of asymptotically conical magnetic surfaces in which current-carrying cylindrical ones are 
nested. When I 00(b) approaches zero, we shall assume that it does so only very slowly, so 
that a WKBJ approach, which considers 1 00(b) as almost constant over large domains of b 
values, will be possible. This approach assumes that the flux distribution on an orthogonal 
trajectory b changes only very slowly when b increases. The same WKBJ treatment can be 
used when -Zoo (6) approaches a non-zero value, since it does so by only slowly varying. As 
explained before, considering this variation is a useful refinement. It is consistent to WKBJ 
analyze Eq.(15) without taking into account any of the other terms present in Eq.(13). The 
gas and poloidal magnetic pressure terms have been shown to be negligible to a certain 
order in r~ l . The relative order of magnitude of the inertia force associated to the curvature 
of the poloidal motion will define the order to which the WKBJ solution is consistent a 
posteriori. It will be shown that this term is indeed insignificant, even when compared to 
the gas pressure term, which is itself small in the field-region. 



4.4. Solution in Field Regions 

Orthogonal trajectories to magnetic surfaces are approximately circles centered at the 
origin. The distribution of flux on such a circle of radius b can be represented by the 
distribution of latitude ip with flux a. In the WKBJ approximation, this distribution slowly 
changes from one circle to the next, so that ip depends not only on a but also weakly on b, 
which we may now identify with the spherical distance R since the orthogonal trajectory's 
label in fact coincides with the distance to the origin all along it. Therefore, magnetic 
surfaces may be locally approximated by cones of semi-opening angle (vr/2 — ip(a,b)). The 
equation for these cones is: 

z — r tan(ip(a, b)) (31) 

From this we obtain by differentiation and ignoring the WKBJ dependance on b: 

1^ , cosib , nnS 

Va = -f— 32 

r dip /da 

From Eq.(25) it results that ip(a) satisfies the differential equation: 

1 dip Vt(a) 



_ _ (33) 

cos^ da /_i I(a, b)^2^E(a) - I(a, b)Q(a)/a(a) 

For this relation to give the flux distribution, the variation with a of I (a, b) must be known, 
which is the case only in field-regions of the asymptotic domain, where I(a, b) becomes 
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independent of a. In a current-carrying boundary layer I(a, b) locally depends strongly on 
a and should be determined by an analysis of the transfield equilibrium. In a field-region, 
where I(a,b) reduces simply to I 00(b), the solution of Eq.(33) is: 

tan(^M)) =tan^(a 1 ,fe))+sinh ( I * \ SM^= ===) (34) 

where a\ is a reference flux in the cell under consideration. If the cell begins at the equator, 
a\ is the flux variable, A, for the equatorial surface and tan(ip(A, b)) = 0. 



4.5. Flux Distribution in Cylindrical Regions of the Field. 

In a region of the free field where the distribution of flux is cylindrical, orthogonal 
trajectories are better represented as straight lines perpendicular to the axis and Eq.(25) 
integrates to: 

roo(a) = rM exp ^£ — - -^========^j (35) 

where a\ is a reference flux in the cylindrical field-region. The flux distribution described by 
Eqs.(34)-(35) with non-zero 1^ is represented, for arbitrarily chosen functions a, E and Q, 
in Fig. (6). A necessary, but not a sufficient, condition for such a solution to be obtained at 
very large distances from the source is that the function (otE/VL) has an absolute minimum. 
Whether a particular system, such as an accretion disk launching a centrifugal wind, can 
indeed meet this necessary condition can only be decided by solving the regularity conditions 
that this flow should satisfy. This point is addressed in the accompanying paper (Heyvaerts 
& Norman 2002b), where we discuss whether the condition that (aE/Q) has a minimum 
value is sufficient to induce cylindrical collimation at infinity. We show that, in a purely 
mathematical sense, it is not but that in a physical sense it should be met in jets of a finite, 
albeit long, extent. 



5. The Polar Boundary Layer 

5.1. Solution in the Polar Boundary Layer. 

The field-region solution does not apply to boundary layers. Within these boundary 
layers, the Lorentz force almost vanishes and forces that would elsewhere be negligible should 
be locally taken into account. The discussion of section (3.1) has shown that the pressure 



-17- 



needs to be considered. In the case of pressureless winds, the poloidal magnetic pressure 
would be the dominant extra force. This situation is considered for completeness in section 
(5.5). Gravity declines to zero at large distances and the azimuthal velocity becomes very 
small at large t/ta- Only if the ratio r/r^ does not become very large would the centrifugal 
force play a role in the transfield equilibrium. In this section, we retain only gas pressure 
and hoop stress in our discussion of transfield equilibrium. 

Physically, the polar boundary layer then locally has the structure of a column pinch. 
The asymptotic transfield equation can be written as Eq.(13), modified by using Eqs.(2) and 
(22), which first gives 



ds a |Va| \ /i a J p |Va 

We find below that the axial density drops down with distance only very slowly, so that the 
poloidal curvature inertia term on the l.h.s. can also be neglected, being smaller than both 
the hoop stress and the pressure term. The transfield equation simplifies to: 

9i Va.V ( ^) + - Va • V (Qp>) = (37) 
a V Po a J P 

Eq.(37) is readily integrated when the region of non- negligible pressure encompasses 
little enough flux that the first integrals E, Q, a and Q can be treated as being constant in 
it, with values Eq, Qq, «o, say. In appendix F, we show this to be so when the asymptotic 
Poynting flux is small compared to the kinetic energy flux. If it is not, the pressure-supported 
region remains similar to a column-pinch, but its structure would now have to be calculated 
numerically. Assuming locally constant first integrals, then, equation (37) becomes 

Va.V (e^ + JL- Q p^) = (38) 
V /"otto 7-1 / 

which integrates, following an orthogonal trajectory to magnetic surfaces, into: 



+ r^T QoP 7 " 1 = t^t QopV (b) (39) 



Po®o 7 - 1 7 

where b is the radius of this quasi-circular orthogonal trajectory and po(b) is the axial density 
at the distance b from the source. Equation (39) can be solved for r in terms of the parameter 

x = p/po (40) 



At the polar axis x — 1 and far from it x decreases to very mall values. Since r = bcosip, 
this solution for r(x) gives tp(x) at a given 6, as expressed by Eq.(43). Close to the polar 
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axis, the Bernoulli equation (5) reduces in the asymptotic domain to: 



V -^=E - -^Qpl 1 (41) 



The Poynting flux is negligible because rB$ vanishes proportionally to r 2 . On the other 
hand, |?7p| is related to |Va| by Eqs.(l) and (2). With the approximation of a constant a(a) 
in the polar boundary layer, this gives, using Eq.(32): 

i-i da/dx 

b 2 p (b)x cos-0 (dip/dx) 

By using Eq.(43) for ip(x), Eqs.(41) and (42), a simple differential equation for a(x) is 
obtained, which is valid for (n/2 — ip) <C I. Its solution is given by Eq.(44). This provides 
the following parametric representation of the solution very near the polar axis: 

^■ ^y d-^) <43) 



a = 



5.2. Matching the Polar Boundary Layer Solution to the Outer Solution. 

Eliminating x, when small (which corresponds to the outer regions of the polar boundary 
layer) between equations (43) and (44) gives the relation between a and ip valid in these 
intermediate regions: 

m ^ = jL-M*r°w exp f a b-p >/2«8 | (45) 



7-1 Q 2 b 2 



7 QoPo 1 (b)p a ^E - ^yQoPo 1 ( b ) , 



It can be asymptotically matched to the outer solution (34). In the vicinity of the polar axis, 
tan^ is large, and the constant term tan^(ai) can be neglected if the reference magnetic 
surface a\ is not one of cylindrical geometry. The relation (34) can be expressed as: 

cos ip = (cosh%) 1 (46) 

where 

X(a, b) = —= / v ' = (47) 

V2p Ioo(b)Ja y/EW-lMQW/aia') 



Near the polar axis x is large and approximately given by 



Q(a')da' 



fin Q> 



VWoo(&) \J ^E{a') - J 00 (6)fi(a')/«(a') y/E Q - I oo (b)Q /a 0i 
which, by Eq.(46), gives the inner limit of the outer solution as: 



(48) 



cos 2 ip — 4 exp — 



V2 



/Voo(fr) 



n(a')da' 



fin a 



(49) 



5.3. Bennet Pinch Relation 

For the exponential arguments in Eqs.(45) and (49) to coincide, it is necessary that: 

7 QoPo (6) = (50) 

7-1 a 

Eq.(50) expresses a relation between the total current supported by the polar boundary layer 
and its inner pressure. The existence of a relation between the total current and the central 
pressure is common to all cylindrical plasma pinches and is usually refered to as a Bennet 
relation. 



5.4. Polar Boundary Layer Current, Density and Radius 

For smooth asymptotic matching, the factors in front of the exponential functions in 
Eqs.(45) and (49) must coincide. Taking Eq.(50) into account, this condition can be written 

as: 

iQo ^alpl- 2 {b) = 4 ^ /_ y/2( 7 - 1) fi n fi(a') da' 

(7-1)^ ~ eXP [ 7/ioa QoPcT 1 (&) Jo V E ( a ') - TJV)W)JW) 

It can be more conveniently written by defining a length £, a dimensionless measure of the 
density, n (b), and a reference magnetic flux a by: 

e = J^)^f- n (b)=po(b)/ PAO a = i/io«ov / 2^^ 2 

(52) 
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and by using the notation 

a = r * (53) 

Jo n y/E(a') - J 00 (6)fi(a')/a(a / ) °o 

The integral A depends on po(b), or n (b), because of Eq.(50). The logarithm of Eq.(51) then 
takes the form: 

Since n (b) is small and (b/£) large, an approximate solution can be obtained by iteration, 
initially neglecting the In (ra ) term on the r.h.s. as compared to n ( - 7 ~ 1 ' ) on the l.h.s.. This 
gives at the simplest degree of approximation: 

^ 1(b) = ^f) (55) 

Eqs.(54) and (55) are dominated by the growth of the logarithm term on their r.h.s. They 
can be satisfied for large 6's in two different ways, according to whether the current Ioo(b) 
at the boundary layer's edge approaches a finite value or decreases to zero. 

When loo (b) approaches a finite value, the Bennet relation (50) shows that the axial 
density should be independent of distance. The logarithmic term in the denominator of 
equation (55) should then be compensated by a divergence of the numerator term A (b) 
(Eq.(53)). This shows that, as b increases, ioo(fr) should approach a limit that causes the 
integral on the r.h.s. of Eq.(53) to diverge. This implies that the asymptotic limit of Ioo{b) 
be the absolute minimum of the function (aE/Q), which, given the presence of a square root 
denominator in the integral on the r.h.s. of Eq.(53), can only be approached from below. 
Therefore Ioo(b) should in this case asymptotically grow towards the absolute minimum of 
(aE/Q). 

If Ioo(b) declines asymptotically to zero, the function A (6) approaches a limit independent 
of b. This indicates that equation (55) should be satisfied by the decrease of n (b) to zero at 
large distances. This is consistent with the fact that in this case all the magnetic surfaces 
flare out parabolically (Heyvaerts & Norman 1989) Equation (55) implies a specific law of 
decrease of po(b), namely, denoting the limit of A (6) for zero current by An: 



An 



7-1 



It is then seen that pn(^) scales with b as (ln(b/£))~~ . The residual current /^(fe), given 
by equation (50), slowly decreases as (ln(26/£)) _1 . 



- 21 - 



The radius r of the circumpolar current channel at distance b is the value of r = b cosijj 
given by Eq.(43) for some intermediate value, of order unity, of the parameter x. This gives 

r ° ~ (7 - 1) n§ (57) 

For Poynting jets, this expression gives the specific value of this radius, since p is given by 
Eq.(50). For kinetic winds, the boundary layer radius slowly increases with distance b as 

(\n(b/£))^h. 

At this point the solution near the pole and in the field region extending in the polar- 
most cell is completely determined. In particular, the solution for the flux distribution in 
the field-region of this cell can, by using Eqs.(49), (51) and (47), be expressed as 

and po and Too are related by Eq.(50). 



5.5. Force- Free Polar Boundary Layers 

If the poloidal magnetic pressure dominates plasma pressure, the corresponding cylin- 
drical structure is described by the transfield mechanical balance equation (13): 



- 2 VL Va ^ / Va 2 
Poa;/ |Va| \2fi r 2 J ' ct |Va 

This equation should be associated with the cylindrical asymptotic form of Eq.(5) 

v 2 P rVLB e , 

E W = ^r- — — , (60) 
2 p: a 

with the mass conservation equation (2), with Eq.(ll) and with Eq.(l). The latter reduces 
in this geometry to 

1 da . . 

Bp = - 1 - 61 

r ar 

The variables vp and p can be eliminated by Eqs.(2) and (11) in favour of the field component 
Bp and the quantity /. Using Eq.(60), B P can then also be expressed in terms of /. We are 
left, for given b, with a pair of ordinary differential equations for I(r) and a(r) where the 
first integrals, which are known functions of a, also appear. This system is 

W ^)+^W , 62) 



dr \r 4 fi 2 \ a J J r 2 dr 
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i *, = v^/ rz™ (63) 

r dr r A \l V ct 
As in the case of pressure-supported polar boundary layers, we assume that the first integrals 
do not vary much over the current carrying boundary layer and treat them as being constants, 
a , Eq and Q say. This is valid for QI/aE <C 1 (app.F) and allows to transform Eq.(62) to 
the simpler form: 

E d fl 2 \ 1 d fl 2 \ n 

^bh^Ur (64) 

This shows that the natural unit length r in the force free pinch, or in other words the core 
radius of this pinch, is: 

< = § < 65 > 

This relation replaces equation (57) for pressure-supported axial boundary layers. The so- 
lution of equation (64) is 

'"■•»> = «FtW (66) 

The solution (66) can be asymptotically matched with the field- region solution (34). Actually 
the solution (66) itself explicitly shows this continuuous transition from I — to I — 1^. It 
establishes the relation, which replaces Eq.(50), between the total electric current carried by 
the polar boundary layer and the poloidal magnetic pressure that supports the associated 
pinching force. The poloidal field and plasma density can be expressed from Eqs.(2) and 
(11) in terms of / as: 

— = ~w - {-') (67) 

and 

p„(6) = ^ Mm { ') (68) 



The limiting value involved is obtained from Eq.(66): 



l im ( £] = 1^1 ( 69 ) 
r-olr 2 / 2Ea y ' 



Equations (67), (68), (69) can be synthesized in the relation: 



B 2 P0 (b) = Ljbp , 
HoPo(b) a 



o 



(70) 



which is the form of Bennet pinch relation appropriate to this case. It differs from Eq.(50) 
only by the substitution of the axial Alfven speed to the axial sound speed. 
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6. Null Surface Boundary Layers 

6.1. Divergence of the Mass Flux to Magnetic Flux Ratio at Null Surfaces 

At a null surface, labelled by a n say, the mass flux to magnetic flux ratio, a (a), defined 
by Eq.(2), diverges as a — > a n if there is mass flux on this null surface. We show in appendix 
(G) that the divergence of a(a) is as 

a(a) ~ ] —— (71) 

\an-a\ v 

where v is a positive number stricly less than unity. Usually, v = 1/2. 

6.2. Structure of Null Surface Boundary Layers. 

The structure of the flow near a null-surface can be derived from the transfield equation 
(13). The ratio of the gas pressure to the toroidal magnetic pressure decreasing to zero as 
the axial distance r increases, the thickness of the pressure-dominated region about the null 
surface, 77, becomes small at large distances. The gradient operator normal to the magnetic 
surface, Va • V/|Va|, noted as — n ■ V, is then of order r\~ x when acting on / or on Qp 1 . 
The left hand side of equation (13) being smaller than (vp/r) (see section (3.1)), can thus 

— » 

be neglected. Similarly the variable r can be treated as a constant because (ft- Vr), of order 
unity, is negligible to (r/77). Eq.(13) then reduces, using Eq.(8), to: 

s '^+S =o (72) 

which, by using Eq.(ll) and regarding Q and Q as constants, Q n and fl n , takes the form 



Since a diverges at a null surface (section 6.1 and app. G) this integrates as 

1 o 2 r 2 Q 2 

QnP 1 + - P —r = QnPlir) (74) 

2 Poa z 

where p n (r) is the density at the cylindrical distance r of the axis on this null surface. Let 
us introduce the parameter 

X = p/p n (r) (75) 
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The mass to flux ratio a (a) is given in terms of X by Eq.(74), namely: 

_J_ = 2QnPr a (r) (±-J-) (76) 

Note that, as it should, a -1 vanishes as the neutral surface is approached, when X approaches 
unity. When the functional dependence of a (a) on a is known, Eq.(76) gives a in terms of X. 
Since this dependence is however not universal (see Eq.(71)), this step cannot be performed 
in a general way. Close to a neutral surface, the Bernoulli equation (5) reduces, in the 
asymptotic domain, to : 

E = f + ^ W - ^ (77) 
2 7 — 1 /i o; 

where can be obtained from Eq.(ll). Using Eq.(76) to account for the variation of a(a) 

in the neighborhood of the neutral surface, we obtain, using Eq.(75) and neglecting the 

variations of r accross the boundary layer, 

V \ » E n - ^—Lq^X^ 1 - 2Q nP l- 1 X- 1 (78) 
2 7 — 1 

On the other hand, \vp\ is related to \Va\ by Eqs.(l) and (2). Using Eq.(32), this provides 
a differential equation for the dependance of the latitude angle ip on a, or equivalently, on 
X, at a given b. Using Eq.(78) this differential equation can be written as: 

, , , a(a)da 

bdijj = v ; (79) 



V2 rX Pn - 2 -^Q nP l' l X^ - 2Q nP r l X-^ 



Eq.(79) could be brought to quadratures for ip(X) if the explicit dependance of a on X could 
be deduced from Eq.(76) for known, and invertible, a(a), such as for example in Eq.(71). 
This step will not be taken here explicitly because the relation between a and a does not have 
an universal character, even near a neutral layer. Although we have not derived the solution 
explicitly, we can still proceed to deduce the necessary conditions for a smooth matching to 
the solution in the far field. 



6.3. Matching the Null Surface Boundary Layer Solution to the Field 

The solution in the boundary layer about the neutral surface is now given by equations 
(76) and (79). It can be asymptotically matched to the field-region solution which is ex- 
pressed in differential form by equation (33). In the field- region near the neutral surface 
Eq.(33) reduces to: 

# = , nnda (80) 
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On the other hand, eliminating X between Eqs.(76) and (79) in the small X limit, which is 
relevant to the outskirts of the equatorial boundary layer, we obtain a differential equation 
for ip valid in this region: 

# = (81) 

ry/2y/2lM)Qnf&(r)yjE n - ^V^/WnrVM 

Matching requires that equations (80) and (81) be identical. This implies that the total 
current at the edge of the null surface boundary layer at a distance b from the wind source 
(corresponding to an axial distance r) is related to the density at the center of the layer by: 

/i /oo(r) = y/2ti Q n r 2 pZ(r) (82) 

This relation expresses the balance between gas pressure at the null surface and the magnetic 
pressure just at the outer edge of its boundary layer, as expected for a sheet-pinch. As a 
result we find that, for Poynting jets, the equatorial density decreases as 

_2 

Pn(r) ~ r -r (83) 

For kinetic winds, ioo(r) decreases as l/(ln (2r /£ cos ip n ))) and the density at the null mag- 
netic surfaces declines as: 



2 



*M~i r H^)r <84) 



6.4. Flux and Current Distribution Near a Neutral Surface. 

Eq.(33) indicates that (dip /da) becomes infinite at a neutral magnetic surface, where I 
vanishes. Does that imply that tan ip becomes infinite at neutral surfaces? In a given current 
cell, Eq.(33) integrates similarly to Eq.(34) to: 

tan(VM)) =tan(^a 1 ,6))+sinh ( H — 1 ,w Ht ) (85) 

\J a V2fx I(a', b) y/E(a') - I(a', b)Q(a')/a(a') J 

Whether tan(ip) diverges or not as a approaches a n depends on the behaviour of the integral 
on the r.h.s. of Eq.(85) as the neutral surface is approached and this in turn depends on how 
I(a, b) varies with (a — a n ) as the neutral surface a n is approached. This can be deduced 
from the solution expressed by Eq.(76). The current I in the asymptotic domain is given by 
Eq.(22). Using Eqs.(75) and (76) it is found that: 

2Q n r 2 pl(r) _ 
Po 
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The parameter X is related to a, or (a — a n ), by Eq.(76). It is shown in section (6.1) that 
a scales with (a n — a) as: 

a' 1 ~ (a n - af (87) 

where v is a positive exponent strictly smaller than unity, usually equal to \. Very near the 
neutral surface, X is close to unity. From Eq.(76), it is found that: 

(fl - a " )2 ^T^ (1 - X7) (88) 
Comparing this with Eq.(86), we conclude that, very near the null surface, 

J~ (l-X r ) 1/2 ~ {a-a n ) v (89) 

With such a dependence of / on a the integral in Eq.(85) converges as a approaches a n , since 
v is less than unity. Therefore, neutral surfaces do not become vertical when approached 
from a conical region. 



7. Shape of the Magnetic Surfaces 

We have now obtained a complete solution in the asymptotic domain, both in field- 
regions (section (4.4)) near the pole (section (5.1)) and near neutral magnetic surfaces (sec- 
tion (6.2)). Since the integral in Eq.(85) converges at a neutral surface, integration may be 
started at the equator, irrespective of whether or not it is a neutral surface. At the equator, 
tan^ vanishes, and thus the integration constant of Eqs.(34) and (85) can be taken to be 
zero. The solution is then extended to other angles ip by using either Eq.(85) or appropriate 
solutions in the neutral or polar boundary layers to specify how J(a, b) depends on a at given 
b. 

The dependence of Ioo{b) on b in field- regions is determined by solving Eq.(54), X(b) 
being defined by Eq.(53) and Iqq by Eq.(50). Eq.(54) may have one or two solutions at large 
b, depending on the function (aE/Q). 

We now have gathered all the information needed to calculate the asymptotic shape of 
magnetic surfaces, both in the field-region and in the various boundary layers. 



7.1. Magnetic Surfaces in Field-regions of Poynting Jets 



When the flow is a Poynting jet, its magnetic structure consists of cylindrical surfaces 
nested into conical ones. The relation between flux and radius for cylindrical surfaces is 



-27- 



given by equations (43) and (44) in the polar boundary layer, noting that r = bcosip, and 
by (58) and (50) outside of it. The dependence of the latitude angle of conical magnetic 
surfaces on flux is given by (34) in the field-region. When there is only one cell extending 
from pole to equator, a\ can be taken as the equatorial flux A and t&nip(a,i) reduces to zero, 
since there is vanishing flux left in the equatorial boundary layer at large distances. The 
shape of the magnetic surfaces in the equatorial boundary layer itself is obtained in section 
7.4. 



7.2. Magnetic Surfaces in Field-regions of Kinetic Winds 

The magnetic surfaces of kinetic winds are described by Eq.(34), Ioo(b) being now given 
by the analysis of section (5.4) (Eqs.(54)-(55)). By Eq.(50) we get, for small Ioo(b): 



7§|J kit (90 > 



We consider the case of polar and equatorial boundary layers of kinetic winds in sections 
7.3 and 7.4 respectively. In the field-regions the shape of magnetic surfaces is given by the 
solution of the following differential equations for r(b) and z(b): 

dr = cos ip (a, b) db dz = sinip (a, b) db (91) 

which are to be integrated in b for constant a. The latter argument will be omitted below. 
The angle i[)(a, b) is given by Eqs.(46)-(47). Since magnetic surfaces are in this case parabolic, 
ip is close to 7r/2 and the spherical distance b can be identified with z. Eq.(91) then reduces 
to 

dr = cosip(z) dz (92) 
Similarly, Eq.(46) for ip simplifies, for large x's and for ioo(&) as given by eq. (90), to: 

/ P \ fc(a) 

cos ^ z)=2 [Yz) (93) 

where k(a) is: 

fA ^(a 7 ) da' 
J a 



rA 

Noting q(a) = 1 — k(a), the solution of Eq.(92) is: 



r 1 (2z^ q{a) 



i q(a) V t- 



(95) 
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The magnetic surfaces are then, in the paraxial region outside of the polar boundary layer, 
a collection of nested power-law paraboloids of variable exponent. Some sections of the 
magnetic surfaces, though extending out of the equatorial boundary layer, may still be close 
enough to the equator that the paraxial approximation if m 7r/2, is inappropriate for them. 
Their shape should be found by integrating Eqs.(91) with no further approximation, as done 
in appendix H. 



7.3. Magnetic Surfaces in the Polar Boundary Layer 

In the polar boundary layer of Poynting jets, magnetic surfaces are cylinders. For kinetic 
winds, the paraxial approximation is fully justified in this region and the shape of magnetic 
surfaces is given by Eqs.(43)-(44), with the notations of Eqs.(52). Eq.(43) can be written as: 

r 2 1/1 1 \ 

(96) 



£ 2 nl'\b) \x x 2 ^ 
whith Eq.(44) providing the value of the parameter x in terms of a by: 

l„I-^(l-^-i) = ^^_ (97) 
x 7-1 V ' ao nl~\b) V ; 

The dimensionless axial density n (b) is approximately given by Eq.(55). The spherical 
distance to the origin being almost identical to z, Eqs.(96)-(97) constitute a set of coupled 
equations relating x, r and z. With the notations of Eqs. (52)- (53) they can be written as: 



r 2 fl 1 \ / In (2z/£) 



2-7 

P -\ x ^A^^r' <98> 

This system can be solved to give r and z in terms of x. For small x we recover the paraxial 
field-region solution, while for x close to unity we obtain the shape of magnetic surfaces in 
the pressure-dominated region very near the polar axis. In this region the magnetic surfaces 
switch from algebraic paraboloids to exponential ones, their shape being given by: 



r fa /2\ 2 ^- 1 ) / 2z\ 2 ^- 1 ) 

~l ~ V On ( A ) k T (100) 
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7.4. Magnetic Surfaces in the Equatorial Boundary Layer 

Let us assume for simplicity that the only neutral surface is the equatorial plane. The 
information on the shape of magnetic surfaces in its boundary layer is provided by the 
parametric solution of Eqs.(76) and (79). The density in the equatorial plane at the distance 
r = b from the source, p e (b), is related to the polar residual current Ioo(b) by Eq.(82). 1 00(b) 
approaches a constant value for Poynting jets and decreases logarithmically for kinetic winds. 
Both cases can be unified by writing 

where m = for Poynting jets and m = 1 for kinetic winds, the factor J m being different in 
the two cases. Equation (82) then gives for the equatorial density 

^={^T {i^kwT (102) 

For small X we recover from Eqs. (76) and (79) the results valid in the field-region. Elimi- 
nating a from Eq.(79) and using Eqs. (76) and (82) we get: 

z = b [ A - (103) 
J a ^ Q e Vpl(b) y/2 (E e - I^nja) 

When loo approaches a non- vanishing constant so does, from Eq.(82), b 2 pj(b). Then z/b 
approaches, at fixed a, a constant value: the magnetic surfaces become conical at the out- 
skirts of the equatorial boundary layer, as they should. For kinetic winds, Eq.(103) gives, 
considering Eq.(lOl): 

tt e (A-a) , , {2b\ ,^ n4S 
z^ v bin — (104) 



J \/2E e 

and the magnetic surfaces become slightly convex paraboloids at the outskirts of the equa- 
torial boundary layer. By contrast, in the region of the equatorial boundary layer where the 
gas pressure dominates, X is close to unity and the integration of Eqs. (76) and (79) gives: 

a(a)da 1 



where again p e (b) is given by Eq.(102). It is then found that z(b) scales, at fixed a, as 

z ~ b^r (106) 
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for Poynting jets and as 




(107) 



for kinetic winds. Note that in this region magnetic surfaces are concave and bend towards 
the equator. This agrees with the conclusions of Okamoto (1999), although the force balance 
considered by this author is between the Lorentz force and the poloidal curvature inertia 
force, whereas we consider balance between the Lorentz force and the pressure gradient 
force. This does not contradict the results of Heyvaerts & Norman (1989) because this is by 
no means the terminal shape of these magnetic surfaces. Indeed, as discussed in section 7.5 
below and in Appendix I, any magnetic field line eventually escapes the equatorial boundary 
layer region, first joining a region at its outskirts where its shape becomes conical or parabolic 
as indicated by Eq.(103) and eventually reaching the field-region. 



7.5. Exit from the Equatorial Boundary Layer 

When magnetic surfaces exit the equatorial boundary layer, they do so with an angle to 
the equatorial plane that decreases with distance. The discussion of section (4.2) assumed 
trajectories orthogonal to magnetic surfaces to cross normal to the equator. For dipolar 
symmetry, an equatorial boundary layer is always present between the field-region and the 
equator itself. The boundary conditions used in section (4.2) are thus consistent only if 
the latitude angle of field lines at the outer edge of the equatorial boundary layer becomes 
increasingly negligible with distance. We show in appendix I that the slope (dz/db) of the 
exiting field line indeed decreases to zero with distance along the equatorial boundary layer. 



7.6. Justification of WKBJ Treatment 

Our WKBJ treatment of the field-region solution is valid only if the inertia force asso- 
ciated with the curvature of the poloidal motion remains negligible. In the case of Poynting 
jets, the poloidal field lines in the field-region asymptotically become exactly straight so that 
solving Eq.(13), supposedly valid to order r^/r, while ignoring the curvature term at the 
l.h.s. of Eq.(12) is obviously consistent. In the case of kinetic winds, the poloidal field lines 
are described by Eqs.(H-5)-(H-6) of appendix H, where their radius of curvature has also 

1 + fc 

been calculated. This radius is proportional to (r~ /k) on field lines for which the exponent 
defined by Eq.(94) is k. Since k approaches unity at the pole, the neglect of the force due 
to poloidal curvature is fully justified in these regions. Near the equator, where k approches 
zero, the curvature radius comes closer to scaling as r, but remains still much larger than r. 
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This is because k never reduces exactly to zero and because of the presence of the dividing 
factor k, which reflects the fact that when field lines become tangent to the equator, their 
curvature must be very small. 

8. Conclusions 

We have derived global solutions for the asymptotic structure of non-relativistic, ro- 
tating, stationary, axisymmetric, polytropic, unconfined, perfect MHD winds. The five la- 
grangian first integrals are assumed to be known. 

The asymptotic structures have been found to consist of vast regions, called field regions, 
which are devoid of any significant residual poloidal electric current density. Residual current 
flows in thin regions in the vicinity of the polar axis and the neighbourhood of null magnetic 
surfaces. Null surfaces can occur at polarity reversals of the wind source or extend over dead 
zones. They delineate separate cells in which the poloidal current achieves closure. 

For kinetic-energy-dominated winds the conversion of total wind energy to kinetic energy 
is shown to progress only logarithmically with distance. 

All winds have been shown to possess a circumpolar current-carrying boundary layer, 
which has the structure of a pressure-supported plasma-jet pinch. Null-surface boundary 
layers have the structure of pressure-supported current sheets. The total electric current is 
constant or slowly diminishes with distance according to an inverse logarithmic law for the 
Poynting flux and kinetic winds respectively. This dimunition is caused by minute amounts 
of current flowing through the diffuse field regions from the pole to the nearest null surface. 

The pressure in the center of these regions, where the toroidal field vanishes, is related 
to the residual current by Bennet pinch relations. The plasma density remains constant at 
the polar axis or declines as a negative power of the logarithm of the distance to the wind 
source as above. Therefore, Poynting flux can be retained, even in kinetic winds, over large 
distances. 

We have calculated the structure of the flow in all possible regions including field regions, 
the polar boundary layer and null-surface boundary layers. The solution is given in terms of 
standard first-integrals using a WKBJ approximation that incorporates the weak dependence 
on the distance from the source. A complete solution has been constructed by asymptotic 
matching of these separate pieces of the solution. Global relations are found between the 
circumpolar current and the density at the polar axis or at neutral surfaces (Eqs.(50) and 
(82)). We have established similar relations in the case of jets with force free polar boundary 
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layers (Eq.(70)). 

The shapes of magnetic surfaces in all parts of the solution and in all relevant regimes 
have been calculated as well. The results are as follows: 

I. For winds which are kinetic-energy dominated at infinity: 

(i) In the free field, the magnetic surfaces focus into algebraic paraboloids (Eq.(95)) as shown 
in Fig. 3. 

(ii) In the polar boundary layer, the magnetic surfaces focus into exponential paraboloids 
(Eq.(lOO)) as shown in Fig.4. 

(iii) Near a null surface, which could be the equatorial plane, the lines are concave, bending 
towards the equator deep inside the neutral boundary layer (Eq.(107)). The magnetic sur- 
faces become straight lines with a logarithmic correction (Eq.(104)) at the edge of the layer. 
These lines are convex, bending away from the equatorial plane, outside the neutral sheet 
as shown in Fig. 5. 

II. For winds carrying Poynting flux at infinity: 

(i) In the free field the solutions asymptote to nested cylindrical and conical magnetic sur- 
faces. (Eqs.(35) and (34)) as shown in Fig.6. 

(ii) In the polar boundary layer they are cylindrical (Eqs.(43)-(44)) as shown in Fig. 7. 

(iii) Near a null surface, which could be the equatorial plane, the lines are concave, bending 
towards the equator deep inside the neutral boundary layer (Eq.(106)). The magnetic sur- 
faces become straight lines (Eq.(103)) at the edge of the layer. These lines remain straight 
outside the neutral sheet as shown in Fig. 8. 

From an observational point of view, the polar and null surface boundary layers which 
carry residual electric current may stand out against the field-regions, both because of their 
large density contrast and because they are a source of free energy. This free energy (as- 
sociated with the currents) has the potential of making them active, by the development 
of instabilities. The density about the pole does not decline with distance in the case of 
Poynting jets, or does it only very slowly in the case of kinetic winds. It may be that what is 
observed as a jet is the dense and active polar boundary layer of a flow developped on a much 
larger angular scale. On these larger scales the flow may be difficult to observe because of 
its very low density and current. Null-surface boundary layers, for example equatorial ones, 
could be observed in association with the jets, although their density and activity is expected 
to decline more rapidly with distance than on the polar axis, because of geometrical effects. 
Their detection could be misinterpreted as accretion disks. 
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The conclusions reached in this paper are of course only valid under our assumptions of 
stationnarity, axisymmetry, polytropic thermodynamics, perfect MHD and unconfinement. 
None of those is expected to be exactly satisfied in reality. Nevertheless we believe that the 
electric circuit picture which emerges from our analysis should be robust against relaxing a 
number of these assumptions. Our description of current closure must remain a feature of 
the solution, because current leakage between main regions of current flow should remain 
weak when the large regions between them suffer little instability. This is expected to 
be the case for regions devoid of current and with only small velocity gradients. Forces 
competing with hoop stress must also remain small under rather general conditions and are 
probably confined to spatially limited regions. By contrast, the distribution of flux in the 
different boundary layers would be modified by the consideration of non-ideal effects, such 
as turbulence resulting from the development of non axisymmetric MHD instabilities. The 
stationnarity assumption is valid on a scale smaller than the size of the cavity carved by the 
wind in the ambient medium. The wind is bordered by a time dependent region featuring 
an expanding shock system. Actually this region is where the residual electric currents flow 
from the circum-polar region to the return regions. That the system, on a scale less than 
that of the global cavity, reaches stationarity seems to be a reasonable assumption. 

In the following paper (Heyvaerts & Norman 2002a) we generalize these results to the 
relativistic regime. An additional paper (Heyvaerts & Norman 2002b) discusses whether 
magnetized outflows are kinetic-energy-dominated or carry Poynting flux. 

The authors thank the Space Telescope Science Institute and the Johns Hopkins Uni- 
versity for continued support to their collaboration. JH also thanks the EC Platon program 
(HPRN-CT-2000-00153) and the Platon collaboration. CN is pleased to thank the Director 
of ESO for support and hospitality during which time this paper was completed. We also 
thank Sundar Srinivasan for significant help with the figures. 



A. Curvature of Poloidal Field Lines 



We write the projection of the equation of motion on the normal to the magnetic surfaces 
as (see for example Heyvaerts & Norman (1989)): 
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where primed quantities are derivatives with respect to a of surface functions. The transfield 
equation (A-l) can be transformed by making use of the following relations which can be 
derived by explicit calculation of (curl v P x v P ) and (curl B P x Bp), using also Eq.(2): 
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The vectors vp and Bp are separated in modulus and direction as vp 
with t the unit vector tangent to the poloidal field line. This gives: 



(A-3) 

vp t and Bp = Bp t, 



(v P ■ V)v P = vp (t- V)f+ t(t- V)(4/2) 



(A-4) 



and a similar equation for Bp. Then use is made of the Fresnet formula (t ■ V)t = n(difj/ds) 
whith ip the angle between the vector t and its projection onto the equatorial plane, s the 
curvilinear abcissa along the poloidal field line and n the unit normal vector oriented towards 

— * 

the polar axis (n = — Va/|Va|). This transforms equations (A-2) and (A-3) into 
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The transfield equation (A-l) is thus reduced to the form: 
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Eliminating vp/2 on its r.h.s. by Eq.(5), Eq.(A-7) becomes: 
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We further eliminate Vg for L and / by using Eqs.(8) and (4) in Eq.(lO), so that 
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This then yields the following form of the transfield equation 
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Eliminating pr 2 for 7 by solving Eq.(A-9), Eq.(A-lO) becomes: 
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The centrifugal term, proportional to da/dr, simplifies by using Eq.(A-9). In the limits 
r — > oo, p/pa — ► 0, $g — > 0, — > (as 1/r), Eq.(A-ll) reduces, by expressing lim(pr 2 ) as 
Poal/Q, to: 

p ds po«/ |Va| / « |Va| V ; 

The terms left in this equation are still not of the same order of magnitude. This is discussed 
in section 3. 



B. Equivalence of Eikonal Equation with Snell's Law 

Let n be the unit vector in the direction of VS, so that the Eikonal equation (28) is: 

VS = N n (B-l) 

The angle of incidence % which appears in Snell's law is the angle between n and VA. Taking 
the curl of equation (B-l) gives Vxn = nx ViV/iV and the projection of this equation onto 
the plane perpendicular to n gives 

(Vxn)xn= — -n — = — (B-2) 

— * 

where the subscript _L indicates the component of VN perpendicular to the local vector n. 
Since n is a unit vector, 

= V(n.n) = 2(n.V)n + 2n x (V x n ) (B-3) 
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This allows to transform the left hand side of Eq.(B-2) into (n.V)n. Introducing the curvi- 
linear abcissa a along an orthogonal trajectory to lines S(r, z) = constant in the meridional 
plane, (n.V)n = dn/da. Thus Eq.(B-2) is: 

* = ^ (B-4) 
da N y ! 

Since n is a unit vector, dn/da is perpendicular to n. According to Eq.(B-4), it is in the 

incidence plane defined by the vectors n and ViV. Let t be the unit vector in the plane of 

incidence which is perpendicular to n and is oriented along the direction of Vj_iV/iV. The 

projection of Eq.(B-4) on the vector t gives: 

t. ^ = = sini |V(lniV)|. (B-5) 
da da 

The minus sign in the second term of (B-5) arises because the curvature of the ray towards 

— * — * 

the direction of ViV causes i to decrease. Since V In N is at an angle i to ft, 

f/lnA |VlniV| cost (I Mi) 



da 

and Eq. (B-5) finally reduces to 



di sini dhvN . 
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da cos i da 
It integrates into Snell's law (ln(iVsini) = constant following a ray). 



C. Field Regions of Winds with a Finite Flux 

If loo is non-zero, the magnetic structure in the asymptotic domain consists of a cylin- 
drically collimated core surrounded by a magnetic structure which encloses poloidal current 
but does not carry it. This structure can be of a conical geometry but it could conceiv- 
ably also be of any current-enclosing parabolic geometry, as for example paraboloids of a 
constant power exponent (Begelman & Li 1994). Our results of section 4.2 show that when 
this magnetic structure is constrained to smoothly join the equatorial plane, the solution 
for flaring magnetic surfaces consists of nested cones. We show here more directly that a 
conical distribution is the only possibility for a current-carrying wind blown by a wind source 
subtending a finite flux. Paraboloids of a constant exponent are represented as: 

z = K(a)r q (C-l) 

with q a constant. The case q = 1 corresponds to cones, while q strictly larger than unity, 
but constant, corresponds to current-enclosing paraboloids. Differentiating Eq.(C-l) gives 

r |Va| = ^\/q 2 + r 2 ( 1 ~i)/K 2 (C-2) 
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Substituting Eq.(C-2) in Eq.(26) with supposedly constant 1^ gives 

K' Q(a) 



Ky/qt + rW-ti/K 2 /ic/ooV^V ' E ( a ) ~ LJMM°) 



(C-3) 



In the large-r limit the l.h.s reduces, when q is strictly larger than unity (paraboloids), to 
a logarithmic derivative. It is therefore impossible in this case to meet the condition that 
K — at the equator, the magnetic surface a = A, if A is finite. This conclusion holds 
also for more complicated current-containing structures of parabolic geometry, such as for 
example z = K(a)e qr , with constant q. In the case of conical asymptotics, we recover the 
results of Heyvaerts & Norman (1989). 



D. Particular Solutions to the Hamilton-Jacobi Equation 

In this appendix we obtain particular, separable, solutions of the Eikonal equation (28) 
r j VjS'I = 1. Squaring it and passing to spherical coordinates, spherical radius R and latitude 
A, it can be written as 

Separable solutions are of the form S(R,X) = F(R) + G(X). The number m being a real 
constant, positive or negative, the variable separation gives for F and G the differential 
equations 

dF 

= m (D-2) 
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d\ J cos 2 A 
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which can be integrated into 



F = m In R + d (D-4) 

g = ± r'dAyi-mwv 

J\ cosA ' 

Simple particular solutions can be found. For m — 0, S depends on A alone, which means 
that flux surfaces which reach r = oo asymptote to cones. For m — ±1, we have two types 
of possible solutions with F = ± In R + C\ and G — ± ln(cos A) + C2, the two ± signs being 
independent of eachother. The only physically acceptable solution is 

S = ± ln(i?cosA) + C (D-6) 
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which gives poloidal field lines (of constant S) parallel to the polar axis. The other sign 
combination must be rejected since it gives rise to magnetic surfaces that do not reach 
z = oo. Equation (D-5) does not integrate simply in the general case. The solution can be 
reduced to quadratures and written as: 

^ , t-, f X dX'y/i — m 2 cos 2 A' _ 

S = C + m InR ± / (D-7) 

Jo cos A' 

For a physical solution, R must be able to approach infinity at constant a, that is at constant 
S. This is possible only when the sign of the functions of R and A on the r.h.s. of Eq.(D-7) 
are different, and the integral in A diverges as R grows very large. For small = n/2 — A, the 
angular integral diverges as ± In 0, and the poloidal field lines approximately become lines 

R \m\ = £05) (D _ 8) 

Since i? ~ z and <fi « (r/z), this gives the cylindrical radius r in terms of z as: 

r = C\S)z {1 -\ m ^ (D-9) 

The physically consistent solutions correspond to \m\ < 1, since otherwise r would not grow 
large with z, contrary to assumptions made to derive this approximate form of the trans- 
field equation. The solution consists in this case of magnetic surfaces of an asymptotically 
parabolic shape, with a power-law exponent independant of the surface. The particular cases 
\m\ — 1 and m = are of course recovered in this more general family of solutions. Note 
however that paraboloids with constant exponent are not valid solutions for a wind source 
subtending a finite flux (Appendix C). Therefore, in this class of separable solutions, only 
those which are asymptotically conical with a cylindrical core of poloidal current-carrying 
magnetic surfaces are acceptable. 



E. Matching across Null Surface Boundary Layers 

We consider here current-carrying jets with several null magnetic surfaces. The function 
Too (6) is non zero, but it has sign jumps at the crossing of the boundary layers about null 
surfaces, which in general are located between pole and equator. The equator itself is, by 
symmetry, a magnetic surface, null or not. Null surfaces separate the space into a number 
of cells. The equatorial plane is embedded in or is bordering one of them. Consider these 
different cells in turn, beginning with one containing the equator, or being contiguous to it. 
From our results of section 4.2, the orthogonal trajectories to magnetic surfaces (rays), are 
circles centered on the polar axis. They must be perpendicular to the equator, which is a 
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particular magnetic surface, and so, they are in this region (asymptotically) centered on the 
origin. The null surface which borders the equatorial cell polewards, being perpendicular to 
these rays, must be a cone centered on the origin. This latter surface is then also perpendic- 
ular to the rays approaching it from the cell polewards, in which 1^ assumes a value which 
differs by its sign from the one in the equatorial cell. The rays in this more poleward cell 
are also circles centered on the polar axis. They must in fact, from this boundary condition, 
be centered at the origin. It is easy to extend this reasoning recurently to show that the 
property of the rays to be centered at the origin passes from the equatorial cell to all the 
other cells. This completes the proof that the magnetic surfaces which reach to infinitely 
large values of r in a Poynting jet must asymptotically be conical. Let us again repeat that 
such a structure must enclose a certain flux about the polar axis in which this finite poloidal 
current is actually flowing, and where poloidal field lines do not follow to infinite values of r. 



F. Flux in the Polar Boundary Layer 

The approximate solution which we have deduced in section 5 has been obtained as- 
suming the first integrals E, Q, Q to be almost constant in the polar boundary layer. The 
validity of this can be judged of by calculating the flux trapped in the polar boundary layer, 
given by equation (44) for, say, x ~ 1/e. This identifies, using Eq.(50), the flux variable at 
the edge of the boundary layer to be: 



a " *> T2 ^7 V - — (F " 1) 

For E, f2, Q to be indeed almost constant in the polar layer, au should be significantly less 
than the total flux A. Note that Ifl/a is of order of the specific energy associated to the 
Poynting flux E p . Then we find 

au . fXpaEpVpo PaE p v qq EpV^ Ep E 

A ~ WA ~ ~ n 2 R 2 A v PA ~ E Q 2 R 2 A v PA [ ' ' 

For fast rotators, the asymptotic velocity is not much larger than vpa, and E is of order 
VL 2 R\. Then the boundary layer flux can only be a small amount of the total flux for flows 
which are largely, but not totally, kinetic-energy-dominated at infinity. 

By contrast, flows which bring a significant Poynting energy flux to infinity would have 
a large part of their asymptotic flux trapped in a high pressure zone of cylindrical geometry. 
Our approximate solution (44)-(43) which assumed little variation of the constants of the 
motion with flux over the polar boundary layer becomes at best schematic in this case. When 
so, the cylindrical pinch equations should be solved exactly in the circumpolar region where 
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pressure is important. So would it also be if (roo/r^) would not be large, as we assumed, in 
a substancial fraction of the polar boundary layer: other forces than gas pressure would then 
have to be taken into account, for example the centrifugal force. The circum-polar structure 
should then be solved for numerically, in the framework of some specific model where the 
first-integrals would be explicitly prescribed. This numerically-determined solution would 
still behave for large r's as Eq.(34). 



G. Mass-to-Magnetic Flux Ratio at a Null Surface 

Near a null magnetic surface of flux variable a = a n the quantity a generally diverges 
as \a n — a| _1//2 , or, more exceptionnally, as some negative power of \a n — a|, with exponent 
smaller in absolute value than unity. This can be seen from Eq.(2). The quantity a, being 
a first-integral, can be evaluated anywhere, for example at a reference sphere of radius Ri 
close to the wind source. If there is to be wind flowing on the surface a n , the mass flux pvp 
should not vanish, while B P vanishes as a — > a n . Let ji{a) be the colatitude at which the 
magnetic surface a is rooted on this reference sphere and let us note the angle /x(a n ) as fj, n . 
Since, close to the source, the magnetic field is essentially undistorted with respect to the 
potential field created by the object at the source of the wind, B P is expected to vanish, 
at given R = Ri, proportionally to \fi(a n ) — fi{a)\ n as fi approaches [i n . The exponent n is 
usually equal to unity and perhaps may sometimes be larger. The flux is the surface integral 
of the component of the field normal to the reference sphere. So we get, for (j,(a) < ji n and 
a <a n , say: 



ft 1 ™ 2ttR 2 sin 

K - o) « / 2nR* sin/i n dp' (/i„ - //)" = ^ ? ^ " M ))"^ (G-l) 

Ja(a) n+1 



For Bp ~ \/i n — /x(a)| n this gives 

Bp w (a n - a)" /(n+1) (G-2) 

Thus a(a) (= pv P /B P ) diverges as (a n — a)~ n ^ n+1 ^> when a — > a n . Since n is expected to 

usually be equal to unity this implies, as expressed by equation (71), that 

a(a)^\(a n -a)\- l l 2 (G-3) 

More generally, with u = (n/n + 1) < 1, 

a{a)^\{a n -a)\ v (G-4) 
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H. Magnetic Surfaces in Field Regions of Kinetic Winds 

The solution in the field-region of kinetic winds is not given accurately enough by (95) 
at non-paraxial surfaces bordering the equator in the field-region. In this appendix we 
calculate their shape which is to be found by integrating equations (91) without any further 
approximation. From Eq.(34), using Eq.(47) and taking the large- x limit we obtain: 

2 

cos ^ — TT~\ ETT (H-l) 

sin ^ = ^j-- {e) ... (H-2) 

where k is defined by Eq.(94). The coordinates r and z are then given in terms of the 
spherical distance b by the differential equations 

dz = (f) fc ~ (f) db 
I ( f )V( f) - i 

— - 2 db 

Eqs.(H-3) (H-4) can be reduced to quadratures to give the following parametric representa- 
tion of these surfaces, taking | = A as the parameter: 



r i r 
1~~k + T Jo 



^ k+1 du 



u (2k/(k+l)) + I ( H " 6 ) 

Again, k is defined by Eq.(94). Let us calculate the radius of curvature of these field lines 
as a function of the parameter A. The line element along them is found to be ds = dX. The 
unit tangent to the line has components 

t - * 2 ^ k m 7) 

r ~ ds~ (2AP + 1 1 } 

tz -Ts~ (2AP + 1 (H - 8) 

and the vector dt/ds is then easily calculated to be equal to n/R c , where n is a unit vector 
perpendicular to t and the curvature radius is 

Rc ~ 4k (2A)*-i (H " 9) 
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Depending on the field line, the parameter k varies between zero near the equator and unity 
at the pole. At large distances r scales as (2A) 1_fe /(l — k) while R c scales as (2A) 1+fc /Ah. 



I. Exit Angle from the Equatorial Boundary Layer 

Let X ex be the value of the parameter X which corresponds to the outer edge of the 
equatorial boundary layer. X ex may be taken as equal to 1/2 or 1/e, say. Assuming a{a) to 
be well approximated in this region by Eq.(71) with v = 1/2, the flux parameter a ex (b) of 
the magnetic surface which exits the boundary layer at b is given by 



^l X exit J \^ J mJ b 



7 



Incidently this shows that the residual flux in the equatorial boundary layer declines with 
distance b. The partial derivative (dz/db) a calculated at a = a ex (b) is the slope of the 
poloidal field line which exits the equatorial boundary layer at b. It is given in terms of b 
and z by 

dz\ 2 — 7 / z N 



(!) <'- 2 > 

\ / a ex 



db J „ 7 

Using Eq.(I-l) and Eqs.(102)-(103), the following expression is obtained: 



dz\ f 2 -l\ 2X2 / 2Q e /i ( 2Q e ((ln(26/£)) m )^ 



S6/o„(6) V 7 J Vm)jW e Xi it \LioJ 2 m J h ^ 



(1-3) 



This slope approaches zero as 6 approaches infinity, as needed for consistency of our analysis 
in section 4.2, in particular for supporting the idea that orthogonal trajectories to magnetic 
surfaces become closer and closer to centered spheres as their radius increases. 



J. Asymptotic Ordering of Lorentz Forces 

All components of B and j approach zero at large distances, though not at the same 
rate. In the field-region, B P) given by Eq.(l), declines with spherical distance as l/b 2 , while 
\B e \ = HqIoq/t may decline slower than 1/6, depending on the shape of magnetic surfaces. 
The toroidal current density is: 

1 -* - I f d Ida d lda\ , T 

3e = — VxB P = + j_i 

Ho Ho \oz r oz or r or j 
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It declines as 1/r 3 , which may be slower than 1/6 3 depending on the geometry of magnetic 

^ — » — * 

surfaces. The poloidal current density jp has a component jp± orthogonal to magnetic 
surfaces and a component jp\\ along them. The latter is constrained by Eq.(15) to vanish, 
keeping the same approximation (or ordering) at which Eq.(15) itself is valid. The discussion 
of section (3.1) has shown that the gradient of gas pressure causes jp\\ to deviate from zero, 
such that 

Q ( aa 
r V b 2 v 



Jp\\Bo « iVQp^l ~ ^ ( — ) (J-2) 



The component jp± is approximately: 

The magnitude of jp±_ depends on how Ioo(b) converges to its limit. If this limit is zero, 
loo is given by Eqs.(50)-(55). If the limit is the minimum I sup of the function (aE/Q), an 
expansion of Eq.(53) for 1^ close to I sup can be made resulting in 

A(6) ~ In [ -7== | (J-4) 

When substituted in Eq.(55) it is found that 1^ is given by 

= hup (^-{iy^j ( j - 5 ) 

where £ is usually a small exponent. It is then deduced that for kinetic winds 

JP± n P(h(f (J " 6) 



while for Poynting jets 

in, ~ J 



Eqs.(J-6) and (J-7) indicate that 

^ ~ p/Tm (J " 8) 

where is a slowly increasing function. We can now describe how the different compo- 

nents of the Lorentz force decline with distance on different magnetic surfaces. The shape of 
these surfaces matters, since the relation between the cylindrical distance r and the radial 
distance b depends on it. It is found that 

jp^Bp ~ ^ (J-9) 
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jp±B e 



1 



(J-10) 



rb*f ± (b) 



1 



(J-ll) 



JeBp 



l 



jp\\Be 



r6 2 T/ 00 (6) 



(J-12) 



The component of the Lorentz force accelerating the fluid along field lines jp±Bg, di- 
minishes most slowly with distance. This is particularly true for kinetic winds, as expected, 
since these transform all of their Poynting energy into kinetic energy. 
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Fig. 1. — Variation of the integrated current in the asymptotic domain following an orthog- 
onal trajectory from the north polar axis to the south polar axis versus the flux variable. 
The poloidal field is assumed to have dipolar symmetry. 
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Fig. 2. — A schematic representation of the poloidal current circuit in the asymptotic do- 
main, in the case of a field with two different polarity zones (dipolar symmetry, left) on the 
wind source (central object), and in the case of three different polarity zones (quadrupolar 
symmetry, right) on the wind source. Current is concentrated in boundary layers near the 
pole and the neutral magnetic surfaces (dotted lines). The heavy arrows indicate the main 
poloidal electric current channels. The dashed line is meant to represent a surface at infinity. 
In the case of kinetic winds there is also a very weak and diffuse current in the field regions 
between the regions of largest current flow so that the boundary layer current vanishes at 
infinity. 
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Fig. 3. — The magnetic field structure for kinetic winds in the asymptotic field region, from 
Eq.(95). The functions E(a) and Q(a) have been taken as constant which implies that 
q(a) = (a/A). On the pole g(0) = and at the equator q(a) = 1. The field lines in each 
quadrant correspond to a/A = 0.2, 0.4, 0.6, 0.8, 0.9 resp. An interpolation formula has been 
used to connect the asymptotic solution to a split-monopole field at the origin. The scale 
for r and z is arbitrary. 
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Fig. 4. — The behaviour of kinetic winds near the polar axis. The panels show: magnetic 
field structure (central frame); the normalized density (Eq.52) at the polar axis vs. the 
distance along the polar axis (upper right); the total current integrated about the polar axis 
vs. the distance along the polar axis (lower right); the ratio of the density to its polar value 
at the same z across the polar pinch (upper left) and the integrated current across the polar 
axis (lower left). The latter two quantities are plotted for some arbitrarily chosen z . The 
central frame uses a value for a such that the proportionality constant in Eq.(100) becomes 
unity. The slightly parabolic shape of the field lines is not clearly visible on the scale of the 
plot. The right panel curves are from Eqs.(50), (52) and (55) with A given by Eq.(53) for 
negligibly small 1^. The left panel curves are from Eqs.(40) (43), (44), (22), (52) and (56). 
The adiabatic index here is 7 = 5/3. 
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Fig. 5. — The behavior of kinetic winds about the equatorial boundary layer. The panels 
show: the magnetic field structure (central); the normalized density n e = p<r '°- > in the equato- 
rial plane vs. distance along the equatorial plane (upper right); the total current integrated 
about the equatorial plane vs. distance along the equatorial plane (lower right); the ratio 
of the density to its equatorial value at the same r across the equatorial sheet (upper left) 
and the integrated current across the equatorial sheet (lower left). The latter two quantities 
are plotted for some arbitrary value of r, r . For simplicity, Q e , a e and Q e have been set 
equal to f2 , a and Q respectively. The field lines (central frame, Eq.(107)) correspond to 
-| = 0.9. Their parabolic shape is not clearly visible on the scale of the plot. The plots in 
the right panels are from Eqs.(82) and (84), while those in the left panels are from Eqs.(75), 
(76), (86) and (87) with v — \. For these plots, the adiabatic index 7 has been taken to be 
equal to |. 
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Fig. 6. — The magnetic field structure for Poynting jets in the field region, from Eqs.(34) 
and (35). From pole to equator, the contours are for (a/A) = 0.4, 0.5, 0.6, 0.7, 0.77, 0.83, 
0.85 and 0.95. The functions E(a) and Q(a) are constants, while ct(a) = (IQ/E){1 + 
[((a/ A) - 0.8) 2 /(l - (a/A)) 1 ' 2 ] }. The unit of the plot for the variables r and z is the scale 
t in Eq.(52). An interpolation formula has been used to connect the asymptotic structure 
of the field to a split monopole field near the origin. Eqs.(34) and (35) become increasingly 
accurate with larger rji and zjl. The scale of the transition between the split monopole 
structure near the origin and the asymptotic solution has been chosen arbitrarily. 
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Fig. 7. — The behaviour of Poynting jets near the polar axis. The plots show: magnetic field 
structure (central frame); the normalized density at the polar axis vs. the distance along the 
polar axis (upper right); the total current integrated about the polar axis vs. the distance 
along the polar axis (lower right); the ratio of the density to its polar value at the same z 
across the polar pinch (upper left) and the integrated current across the polar axis (lower 
left). The latter two quantities are plotted for some arbitrarily chosen z . The central frame 
shows that the field lines are exact cylinders (Eq.(57)). The right panel curves use Eqs.(50) 
and (52). The left panel curves are obtained from Eqs.(40), (43), (44), (22) and (52). /„ 
here is the absolute minimum of 
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Fig. 8. — The behaviour of Poynting jets about the equatorial boundary layer. The plots 
show: magnetic field structure (central frame, from Eq.(106)); the normalized density n e = 
in the equatorial plane vs. distance along the equatorial plane (upper right); the total 
current integrated about the equatorial plane vs. distance along the equatorial plane (lower 
right); the ratio of the density to its equatorial value at the same r across the equatorial 
sheet (upper left) and the integrated current across the equatorial sheet (lower left). The 
latter two quantities are plotted for some arbitrary value of r, r . The field structure (central 
frame) is plotted for ^ = 0.9. The plots in the right panels are from Eqs.(82) and (83), while 
those in the left panels are from Eqs.(75), (76), (86) and (87) with v — \. For these plots, 
the adiabatic index 7 has been taken to be equal to |. 



